#include"adc2.h"
void ADC2::setup2p1h()
{
  states_2p1h.resize(numGroups);
  for(int p1=FermiSurface+1;p1<totalOrbitals;p1++)
    {
      for(int p2=p1;p2<totalOrbitals;p2++)
	{
	  for(int h1=0;h1<=FermiSurface;h1++)
	    {
	      for(int nb=0;nb<numGroups;nb++)
		{
		  int l=blocks[nb].l;
		  int jj=blocks[nb].jj;
		  int tz=blocks[nb].tz;

		  if((HFOrbs[p1].l + HFOrbs[p2].l-HFOrbs[h1].l-l)%2) continue;
		  if( (HFOrbs[p1].jj+HFOrbs[p2].jj+HFOrbs[h1].jj) <jj ) continue;
		  if( abs(abs(HFOrbs[p1].jj - HFOrbs[p2].jj) - HFOrbs[h1].jj) > jj ) continue;
		  if((HFOrbs[p1].tz + HFOrbs[p2].tz-HFOrbs[h1].tz-tz)!=0) continue;
		  int jmax=(HFOrbs[p1].jj+HFOrbs[p2].jj)/2;
		  int jmin=abs(HFOrbs[p1].jj-HFOrbs[p2].jj)/2;
		  for(int J=jmin;J<=jmax;J++)
		    {
		      if( (p1==p2)&& (J%2) ) continue;
		      states_2p1h[nb].push_back(state_2p1h(p1,p2,h1,J));
		    }
		}
	    }
	}
    }
}
void ADC2::setup2h1p()
{
  states_2h1p.resize(numGroups);
  for(int h1=0;h1<=FermiSurface;h1++)
    {
      for(int h2=h1;h2<=FermiSurface;h2++)
	{
	  for(int p1=FermiSurface+1;p1<totalOrbitals;p1++)
	    {
	      for(int nb=0;nb<numGroups;nb++)
		{
		  int l=blocks[nb].l;
		  int jj=blocks[nb].jj;
		  int tz=blocks[nb].tz;
		  if((HFOrbs[h1].l + HFOrbs[h2].l-HFOrbs[p1].l-l)%2) continue;
		  if(HFOrbs[h1].jj+HFOrbs[h2].jj+HFOrbs[p1].jj<jj) continue;
		  if(abs(abs(HFOrbs[h1].jj - HFOrbs[h2].jj) - HFOrbs[p1].jj) > jj ) continue;
		  if((HFOrbs[h1].tz + HFOrbs[h2].tz-HFOrbs[p1].tz-tz)!=0) continue;
		  int jmax=(HFOrbs[h1].jj+HFOrbs[h2].jj)/2;
		  int jmin=abs(HFOrbs[h1].jj-HFOrbs[h2].jj)/2;
		  for(int J=jmin;J<=jmax;J++)
		    {
		      if( (h1==h2)&& (J%2) ) continue;
		      states_2h1p[nb].push_back(state_2h1p(h1,h2,p1,J));
		    }
		}
	    }
	}
    }
}
void ADC2::setup()
{
  totalOrbitals=HFOrbs.size();
  FermiSurface=pSystem_Tz_HF->FermiSurface;
  numGroups=pSystem_Tz_HF->Groups.size();
  blocks.resize(numGroups);
  for(int nb=0;nb<numGroups;nb++)
    {
      blocks[nb].dim=Groups[nb].size();
      const HFOrbital & Orb=HFOrbs[Groups[nb][0]];
      blocks[nb].l=Orb.l;
      blocks[nb].jj=Orb.jj;
      blocks[nb].tz=Orb.tz;
    }
  Z.resize(numGroups);
  poles.resize(numGroups);
  sep_energy=(HFOrbs[FermiSurface].e+HFOrbs[FermiSurface+1].e)/2;
  setup2p1h();
  setup2h1p();
}

void ADC2::solveDysonMatrix()
{
#pragma omp parallel for
  for(int nb=0;nb<numGroups;nb++)
    {
      int dim_static=blocks[nb].dim;
      int dim_2p1h=states_2p1h[nb].size();
      int dim_2h1p=states_2h1p[nb].size();
      int dim=dim_static+dim_2p1h+dim_2h1p;
      MatrixXd mat=MatrixXd::Zero(dim,dim);

      for(int i=0;i<dim_static;i++)
	{
	  int a= Groups[nb][i];
	  mat(i,i)=HFOrbs[a].e;
	  for(int j=0;j<dim_2p1h;j++)
	    {
	      int p1=states_2p1h[nb][j].p1;
	      int p2=states_2p1h[nb][j].p2;
	      int h1=states_2p1h[nb][j].h1;
	      int J=states_2p1h[nb][j].J;
      	      double factor = ( (p1==p2)?sqrt(0.5):1.0 );
	      mat(i,dim_static+j)=sqrt(2*J+1)/sqrt(blocks[nb].jj+1)*pSystem_Tz_HF->get2B(a,h1,p1,p2,J)*factor;
	      mat(dim_static+j,i)=mat(i,dim_static+j);
	      mat(dim_static+j,dim_static+j)=HFOrbs[p1].e+HFOrbs[p2].e-HFOrbs[h1].e;
	    }
	  for(int j=0;j<dim_2h1p;j++)
	    {
	      int h1=states_2h1p[nb][j].h1;
	      int h2=states_2h1p[nb][j].h2;
	      int p1=states_2h1p[nb][j].p1;
	      int J=states_2h1p[nb][j].J;
	      double factor = ( (h1==h2)?sqrt(0.5):1.0 );
	      mat(i,dim_static+dim_2p1h+j)=sqrt(2*J+1)/sqrt(blocks[nb].jj+1)*pSystem_Tz_HF->get2B(a,p1,h1,h2,J)*factor;
	      mat(dim_static+dim_2p1h+j,i)=mat(i,dim_static+dim_2p1h+j);
	      mat(dim_static+dim_2p1h+j,dim_static+dim_2p1h+j)=HFOrbs[h1].e+HFOrbs[h2].e-HFOrbs[p1].e;
	    }
	}      
      SelfAdjointEigenSolver<MatrixXd> es(mat);
      Z[nb]=es.eigenvectors().topLeftCorner(dim_static,dim);
      poles[nb]=es.eigenvalues();
    }
}

void ADC2::solveDysonMatrix_Lanczos(const int p,const int q)
{
  omp_set_num_threads(5);
#pragma omp parallel for
  for(int nb=0;nb<numGroups;nb++)
    {
      int dim_static=blocks[nb].dim;
      int dim_2p1h=states_2p1h[nb].size();
      int dim_2h1p=states_2h1p[nb].size();

      vector<int> n_iter_fw(dim_static,p);
      vector<int> n_iter_bw(dim_static,q);

      int P=p,Q=q;
      int dim_fw=P*dim_static;
      int dim_bw=Q*dim_static;
      if(dim_fw>dim_2p1h)
      	{
      	  P=dim_2p1h/dim_static;
      	  for(int i=0;i<dim_static;i++)
      	    {
      	      n_iter_fw[i]=P;
      	    }
	  n_iter_fw[0]+=dim_2p1h%dim_static;
      	  dim_fw=dim_2p1h;
      	}
      if(dim_bw>dim_2h1p)
      	{
	  Q=dim_2h1p/dim_static;
      	  for(int i=0;i<dim_static;i++)
      	    {
      	      n_iter_bw[i]=Q;
      	    }
   	  n_iter_bw[0]+=dim_2h1p%dim_static;
      	  dim_bw=dim_2h1p;
      	}
      
      int dim=dim_static+dim_fw+dim_bw;
      MatrixXd mat=MatrixXd::Zero(dim,dim);

      MatrixXd M_A_fw(dim_2p1h,dim_static);
      MatrixXd M_C_fw = MatrixXd::Zero(dim_2p1h,dim_2p1h);

      MatrixXd M_A_bw(dim_2h1p,dim_static);
      MatrixXd M_C_bw = MatrixXd::Zero(dim_2h1p,dim_2h1p);

      for(int i=0;i<dim_static;i++)
	{
	  int a=Groups[nb][i];
	  mat(i,i)=HFOrbs[a].e;
	  for(int j=0;j<dim_2p1h;j++)
	    {
	      int p1=states_2p1h[nb][j].p1;
	      int p2=states_2p1h[nb][j].p2;
	      int h1=states_2p1h[nb][j].h1;
	      int J=states_2p1h[nb][j].J;
      	      double factor = ( (p1==p2)?sqrt(0.5):1.0 );
	      M_A_fw(j,i)=sqrt(2*J+1)/sqrt(blocks[nb].jj+1)*pSystem_Tz_HF->get2B(a,h1,p1,p2,J)*factor;
	      M_C_fw(j,j)=HFOrbs[p1].e+HFOrbs[p2].e-HFOrbs[h1].e;
	    }
	  for(int j=0;j<dim_2h1p;j++)
	    {
	      int h1=states_2h1p[nb][j].h1;
	      int h2=states_2h1p[nb][j].h2;
	      int p1=states_2h1p[nb][j].p1;
	      int J=states_2h1p[nb][j].J;
	      double factor = ( (h1==h2)?sqrt(0.5):1.0 );
	      M_A_bw(j,i)=sqrt(2*J+1)/sqrt(blocks[nb].jj+1)*pSystem_Tz_HF->get2B(a,p1,h1,h2,J)*factor;
	      M_C_bw(j,j)=HFOrbs[h1].e+HFOrbs[h2].e-HFOrbs[p1].e;
	    }
	}
      buildLanczosMatrix( mat.block(dim_static,0,dim_fw,dim_static),mat.block(dim_static,dim_static,dim_fw,dim_fw), M_A_fw,M_C_fw, n_iter_fw);

      //it's commented below because we only use the lower triangular part of symmetry matrix "mat" to diagnolize. It's not necessary to give all the elements of "mat"
      //      mat.block(0,dim_static,dim_static,dim_fw)=mat.block(dim_static,0,dim_fw,dim_static).transpose();

      buildLanczosMatrix( mat.block(dim_static+dim_fw,0,dim_bw,dim_static),mat.block(dim_static+dim_fw,dim_static+dim_fw,dim_bw,dim_bw), M_A_bw,M_C_bw ,n_iter_bw);
      //it's commented below because we only use the lower triangular part of symmetry matrix "mat" to diagnolize. It's not necessary to give all the elements of "mat"
      //      mat.block(0,dim_static+dim_fw,dim_static,dim_bw)=mat.block(dim_static+dim_fw,0,dim_bw,dim_static).transpose();

      SelfAdjointEigenSolver<MatrixXd> es(mat);
      Z[nb]=es.eigenvectors().topLeftCorner(dim_static,dim);
      poles[nb]=es.eigenvalues();
      
    }
}

/*
     /          |         \            /          |         \ 
     | M_static |   M_A'  |	       | M_static |     A'  |  
  Q' |__________|_________|  Q =====>  |__________|_________|  
     |          |         |	       |          |         |  
     |    M_A   |   M_C   |	       |    A     |     C   |  
     \          |         /	       \          |         /  

         /          |           \ 
         |     I    |    0      |
  Q  =   |__________|___________|
         |          |           |
         |     0    | Lanczos_V |
         \          |           /

*/
void ADC2::buildLanczosMatrix(Ref<MatrixXd> A,Ref<MatrixXd> C,const MatrixXd & M_A,const MatrixXd & M_C,const vector<int>& n_iter)
{
  int dim_static=M_A.cols();
  int dim_before=M_A.rows();
  int dim_after=A.rows();
  if(dim_after==0) return;
  MatrixXd Lanczos_V=MatrixXd::Zero(dim_before,dim_after);
  MatrixXd U(dim_before,dim_static);
  int index=0;
  for(int i=0;i<dim_static;i++)
    {
      VectorXd V0=M_A.col(i);
      //orthonormalize V0 with Lanczos_V.col( 0,1,2... index-1 )
      for(int j=index-1;j>=0;j--)
	{
	  V0-=V0.dot( Lanczos_V.col(j) )*Lanczos_V.col(j);
	}
      V0.normalize();
      Lanczos_V.col(index)=V0;
      VectorXd V_next=M_C*Lanczos_V.col(index);
      C(index,index)=Lanczos_V.col(index).dot(V_next);
      for(int j=1;j< n_iter[i] ;j++)
	{
	  //ortholize V_next with Lanczos_V.col(0...index+j-1)
	  for(int k=index+j-1;k>=0;k--)
	    {
	      V_next-=V_next.dot( Lanczos_V.col(k) )*Lanczos_V.col(k);
	    }
	  C(index+j-1,index+j)=V_next.norm();
	  C(index+j,index+j-1)=C(index+j-1,index+j);
	  V_next/=C(index+j-1,index+j);
	  Lanczos_V.col(index+j)=V_next;
	  V_next=M_C*Lanczos_V.col(index+j);
	  C(index+j,index+j)=Lanczos_V.col(index+j).dot(V_next);
	}
      U.col(i)=V_next;
      //because we use block-like Lanczos algorithm, the matrix is not tridiagonal, it's fishbone-like, we calculate the matrix element of the wings below. 
      int index_j=0;
      for(int j=0;j<i;j++)
	{
	  for(int k=index;k<index+n_iter[i];k++)
	    {
	      C( k,index_j+n_iter[j] - 1 ) =  Lanczos_V.col(k).dot( U.col(j) );
	      C( index_j+n_iter[j] - 1, k) = C( k,index_j+n_iter[j] - 1 );
	    }
	  index_j+=n_iter[j];
	}
      index+=n_iter[i];
    }
  A=Lanczos_V.transpose()*M_A;
}

void ADC2::calEnergy()
{
  Energy=0;
  num=0;
  for(int nb=0;nb<numGroups;nb++)
    {
      int dim=blocks[nb].dim;
      for(int i=0;i<poles[nb].rows();i++)
	{
	  if(poles[nb](i)<sep_energy)
	    {
	      for(int bra=0;bra<dim;bra++)
		{
		  Energy+=poles[nb](i)*Z[nb](bra,i)*Z[nb](bra,i)*(blocks[nb].jj+1.);
		  num+=Z[nb](bra,i)*Z[nb](bra,i)*(blocks[nb].jj+1.);
		  for(int ket=0;ket<dim;ket++)
		    {
		      Energy+=get1B(nb,bra,ket)*Z[nb](bra,i)*Z[nb](ket,i)*(blocks[nb].jj+1.);
		    }
		}
	    }
	}
    }
  Energy/=2;
}




void ADC2::solveDysonMatrix_Lanczos_iter(const int p,const int q)
{
  vector< MatrixXd > Z_new(numGroups);
  vector< VectorXd > poles_new(numGroups);
  omp_set_num_threads(5);
#pragma omp parallel for
  for(int nb=0;nb<numGroups;nb++)
    {
      int dim_fw_all=0;
      int dim_bw_all=0;
      int l=blocks[nb].l;
      int jj=blocks[nb].jj;
      int tz=blocks[nb].tz;
      for(int nb1=0;nb1<numGroups;nb1++)
	{
	  for(int nb2=nb1;nb2<numGroups;nb2++)
	    {
	      for(int nb3=0;nb3<numGroups;nb3++)
		{
		  int l1=blocks[nb1].l;
		  int jj1=blocks[nb1].jj;
		  int tz1=blocks[nb1].tz;

		  int l2=blocks[nb2].l;
		  int jj2=blocks[nb2].jj;
		  int tz2=blocks[nb2].tz;

		  int l3=blocks[nb3].l;
		  int jj3=blocks[nb3].jj;
		  int tz3=blocks[nb3].tz;

		  if((l1 + l2-l3-l)%2) continue;
		  if( (jj1+jj2+jj3) <jj ) continue;
		  if( abs(abs(jj1 - jj2) - jj3) > jj ) continue;
		  if((tz1 + tz2-tz3-tz)!=0) continue;

		  int dim1=poles[nb1].size();
		  int dim2=poles[nb2].size();
		  int dim3=poles[nb3].size();
		  for(int index1=0;index1<dim1;index1++)
		    for(int index2=((nb1==nb2)?index1:0);index2<dim2;index2++)
		      for(int index3=0;index3<dim3;index3++)
			{
			  int jmax=(jj1+jj2)/2;
			  int jmin=abs(jj1-jj2)/2;
			  if(poles[nb1](index1)>sep_energy&&poles[nb2](index2)>sep_energy&&poles[nb3](index3)<sep_energy)
			    {
			      for(int J=jmin;J<=jmax;J++)
				{
				  dim_fw_all++;
				}
			    }
			  else if(poles[nb1](index1)<sep_energy&&poles[nb2](index2)<sep_energy&&poles[nb3](index3)>sep_energy)
			    {
			      for(int J=jmin;J<=jmax;J++)
				{
				  dim_bw_all++;
				}
			    }
			}
		}
	    }
	}

      int dim_static=blocks[nb].dim;
      vector<int> n_iter_fw(dim_static,p);
      vector<int> n_iter_bw(dim_static,q);
      int P=p,Q=q;
      int dim_fw=P*dim_static;
      int dim_bw=Q*dim_static;
      if(dim_fw>dim_fw_all)
	{
	  P=dim_fw_all/dim_static;
	  for(int i=0;i<dim_static;i++)
	    {
	      n_iter_fw[i]=P;
	    }
  	  n_iter_fw[0]+=dim_fw_all%dim_static;
	  dim_fw=dim_fw_all;
	}
      if(dim_bw>dim_bw_all)
	{
  	  Q=dim_bw_all/dim_static;
	  for(int i=0;i<dim_static;i++)
	    {
	      n_iter_bw[i]=Q;
	    }
     	  n_iter_bw[0]+=dim_bw_all%dim_static;
	  dim_bw=dim_bw_all;
	}
      int dim=dim_static+dim_fw+dim_bw;
      MatrixXd mat=MatrixXd::Zero(dim,dim);

      MatrixXd M_A_fw = MatrixXd::Zero(dim_fw_all,dim_static);
      VectorXd M_C_fw = VectorXd::Zero(dim_fw_all);

      MatrixXd M_A_bw = MatrixXd::Zero(dim_bw_all,dim_static);
      VectorXd M_C_bw = VectorXd::Zero(dim_bw_all);

      int n_fw=0,n_bw=0;
      for(int nb1=0;nb1<numGroups;nb1++)
	{
	  for(int nb2=nb1;nb2<numGroups;nb2++)
	    {
	      for(int nb3=0;nb3<numGroups;nb3++)
		{
		  int l1=blocks[nb1].l;
		  int jj1=blocks[nb1].jj;
		  int tz1=blocks[nb1].tz;

		  int l2=blocks[nb2].l;
		  int jj2=blocks[nb2].jj;
		  int tz2=blocks[nb2].tz;

		  int l3=blocks[nb3].l;
		  int jj3=blocks[nb3].jj;
		  int tz3=blocks[nb3].tz;

		  if((l1 + l2-l3-l)%2) continue;
		  if( (jj1+jj2+jj3) <jj ) continue;
		  if( abs(abs(jj1 - jj2) - jj3) > jj ) continue;
		  if((tz1 + tz2-tz3-tz)!=0) continue;

		  int dim1=poles[nb1].size();
		  int dim2=poles[nb2].size();
		  int dim3=poles[nb3].size();
		  for(int index1=0;index1<dim1;index1++)
		    for(int index2=((nb1==nb2)?index1:0);index2<dim2;index2++)
		      for(int index3=0;index3<dim3;index3++)
			{
			  int jmax=(jj1+jj2)/2;
			  int jmin=abs(jj1-jj2)/2;
			  if(poles[nb1](index1)>sep_energy&&poles[nb2](index2)>sep_energy&&poles[nb3](index3)<sep_energy)
			    {
			      for(int J=jmin;J<=jmax;J++)
				{
				  for(int i=0;i<dim_static;i++)
				    {
				      double temp=0;
				      for(int ia=0;ia<blocks[nb1].dim;ia++)
					for(int ib=0;ib<blocks[nb2].dim;ib++)
					  for(int ic=0;ic<blocks[nb3].dim;ic++)
					    {
					      temp+=Z[nb1](ia,index1)*Z[nb2](ib,index2)*Z[nb3](ic,index3)*get2B(nb,i,nb3,ic,nb1,ia,nb2,ib,J);
					    }
				      double factor=( ((nb1==nb2)&&(index1==index2))?sqrt(0.5):1.0 );
				      M_A_fw(n_fw,i)=temp*factor*sqrt(2*J+1.)/sqrt(jj+1.);
				    }
				  M_C_fw(n_fw)=poles[nb1](index1)+poles[nb2](index2)-poles[nb3](index3);
				  n_fw++; 
				}
			    }
			  else if(poles[nb1](index1)<sep_energy&&poles[nb2](index2)<sep_energy&&poles[nb3](index3)>sep_energy)
			    {
			      for(int J=jmin;J<=jmax;J++)
				{
				  for(int i=0;i<dim_static;i++)
				    {
				      double temp=0;
				      for(int ia=0;ia<blocks[nb1].dim;ia++)
					for(int ib=0;ib<blocks[nb2].dim;ib++)
					  for(int ic=0;ic<blocks[nb3].dim;ic++)
					    {
					      temp+=Z[nb1](ia,index1)*Z[nb2](ib,index2)*Z[nb3](ic,index3)*get2B(nb,i,nb3,ic,nb1,ia,nb2,ib,J);
					    }
				      double factor=( ((nb1==nb2)&&(index1==index2))?sqrt(0.5):1.0 );
				      M_A_bw(n_bw,i)=temp*factor*sqrt(2*J+1.)/sqrt(jj+1.);
				    }
				  M_C_bw(n_bw)=poles[nb1](index1)+poles[nb2](index2)-poles[nb3](index3);
				  n_bw++;
				}
			    }
			}
		}
	    }
	}

      
      for(int i=0;i<dim_static;i++)
	{
	  for(int j=i;j<dim_static;j++)
	    {
	      mat(j,i)=get1B(nb,j,i);
	      double temp=0;
	      for(int nb_temp=0;nb_temp<numGroups;nb_temp++)
		{
		  int jmin=abs(blocks[nb_temp].jj - jj)/2;
		  int jmax=(blocks[nb_temp].jj + jj)/2;
		  
		  int dim_temp=blocks[nb_temp].dim;
		  for(int k=0;k<dim_temp;k++)
		    {
		      for(int l=0;l<dim_temp;l++)
			{
			  for(int J=jmin;J<=jmax;J++)
			    {
			      for(int a=0;a<poles[nb_temp].size();a++)
				{
				  if(poles[nb_temp](a)>sep_energy) continue;
				  temp+=(2*J+1)*get2B(nb,j,nb_temp,k,nb,i,nb_temp,l,J)*Z[nb_temp](k,a)*Z[nb_temp](l,a);
				}
			    }
			}
		    }
		}
	      mat(j,i)+=temp/(jj+1);
	    }
	}
     
      buildLanczosMatrix_iter( mat.block(dim_static,0,dim_fw,dim_static),mat.block(dim_static,dim_static,dim_fw,dim_fw), M_A_fw,M_C_fw, n_iter_fw);

      //it's commented below because we only use the lower triangular part of symmetry matrix "mat" to diagnolize. It's not necessary to give all the elements of "mat"
      //      mat.block(0,dim_static,dim_static,dim_fw)=mat.block(dim_static,0,dim_fw,dim_static).transpose();

      buildLanczosMatrix_iter( mat.block(dim_static+dim_fw,0,dim_bw,dim_static),mat.block(dim_static+dim_fw,dim_static+dim_fw,dim_bw,dim_bw), M_A_bw,M_C_bw ,n_iter_bw);
      //it's commented below because we only use the lower triangular part of symmetry matrix "mat" to diagnolize. It's not necessary to give all the elements of "mat"
      //      mat.block(0,dim_static+dim_fw,dim_static,dim_bw)=mat.block(dim_static+dim_fw,0,dim_bw,dim_static).transpose();

      SelfAdjointEigenSolver<MatrixXd> es(mat);
      Z_new[nb]=es.eigenvectors().topLeftCorner(dim_static,dim);
      poles_new[nb]=es.eigenvalues();
    }
  Z=Z_new;
  poles=poles_new;
}



void ADC2::buildLanczosMatrix_iter(Ref<MatrixXd> A,Ref<MatrixXd> C,const MatrixXd & M_A,const VectorXd & M_C,const vector<int>& n_iter)
{
  int dim_static=M_A.cols();
  int dim_before=M_A.rows();
  int dim_after=A.rows();
  if(dim_after==0) return;
  MatrixXd Lanczos_V=MatrixXd::Zero(dim_before,dim_after);
  MatrixXd U(dim_before,dim_static);
  int index=0;
  for(int i=0;i<dim_static;i++)
    {
      VectorXd V0=M_A.col(i);
      //orthonormalize V0 with Lanczos_V.col( 0,1,2... index-1 )
      for(int j=index-1;j>=0;j--)
	{
	  V0-=V0.dot( Lanczos_V.col(j) )*Lanczos_V.col(j);
	}
      V0.normalize();
      Lanczos_V.col(index)=V0;
      VectorXd V_next=M_C.cwiseProduct( Lanczos_V.col(index) );
      C(index,index)=Lanczos_V.col(index).dot(V_next);
      for(int j=1;j< n_iter[i] ;j++)
	{
	  //ortholize V_next with Lanczos_V.col(0...index+j-1)
	  for(int k=index+j-1;k>=0;k--)
	    {
	      V_next-=V_next.dot( Lanczos_V.col(k) )*Lanczos_V.col(k);
	    }
	  C(index+j-1,index+j)=V_next.norm();
	  C(index+j,index+j-1)=C(index+j-1,index+j);
	  V_next/=C(index+j-1,index+j);
	  Lanczos_V.col(index+j)=V_next;
	  V_next=M_C.cwiseProduct( Lanczos_V.col(index+j) );
	  C(index+j,index+j)=Lanczos_V.col(index+j).dot(V_next);
	}
      U.col(i)=V_next;
      //because we use block-like Lanczos algorithm, the matrix is not tridiagonal, it's fishbone-like, we calculate the matrix element of the wings below. 
      int index_j=0;
      for(int j=0;j<i;j++)
	{
	  for(int k=index;k<index+n_iter[i];k++)
	    {
	      C( k,index_j+n_iter[j] - 1 ) =  Lanczos_V.col(k).dot( U.col(j) );
	      C( index_j+n_iter[j] - 1, k) = C( k,index_j+n_iter[j] - 1 );
	    }
	  index_j+=n_iter[j];
	}
      index+=n_iter[i];
    }
  A=Lanczos_V.transpose()*M_A;
}

void ADC2::printResults() const
{
  cout<<"sep:\t"<<sep_energy<<endl;
  for(int nb=0;nb<blocks.size();nb++)
    {
      cout<<"--------\n";
      cout<<blocks[nb].l<<"\t"<<blocks[nb].jj<<"\t"<<blocks[nb].tz<<endl;
      //      cout<<poles[nb].transpose()<<endl;
      //      cout<<"********\n";
      for(int i=0;i<poles[nb].rows();i++)
	{
	  double norm=Z[nb].col(i).norm();
	  if(norm>0.5)
	    cout<<poles[nb](i)<<"\t"<<norm<<endl;
	}
      cout<<"--------\n";      
    }
}
